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Abstract. Following Toomre & Kalnajs (1991), local models of slightly dissipative self-gravitating disks show how 
inhomogeneous structures can be maintained over several galaxy rotations. Their basic physical ingredients are 
self-gravity, dissipation and differential rotation. In order to explore the structures resulting from these processes 
on the kpc scale, local simulation of self-gravitating disks are performed in this paper in 2D as well as in 3D. 
The third dimension becomes a priori important as soon as matter clumping causes a tight coupling of the 3D 
equations of motion. The physically simple and general framework of the model permits to make conclusions 
beyond the here considered scales. A time dependent affine coordinate system is used, allowing to calculate the 
gravitational forces via a particle-mesh FFT-method, increasing the performance with respect to previous direct 
force calculations. 

Persistent patterns, formed by transient structures, whose intensity and morphological characteristic depend on the 
dissipation rate are obtained and described. Some of our simulations reveal first signs of mass-size and velocity 
dispersion-size power-law relations, but a clear scale invariant behavior will require more powerful computer 
techniques. 
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1. Introduction 

Classical gravity is scale free, i.e., self-gravitating systems 
may form similar structures on different scales. Indeed, ob- 
servations of the interstellar medium, spiral disks and cos- 
mic structures, do reveal similar characteristics. Although 
the structures in these systems are lumpy and inhomoge- 
neous, they do not seem yet completely random. 

The observations of molecular clouds reveal hierarchi- 
cal structures with power-law behavior over several orders 
of magnitude in scale (Larson [1981 , Scalo 1985, Falgarone 
et al. |1992|, Heithausen et al. 1999). Larson (|197"8|) found 



and the cosmic structure is that the matter distribution 
can be characterized by a comparable fractal dimension. 
The cosmic and the interstellar fractal dimensions are, 
-P caiaxic s ~ 2 ± 0.2 (Sylos Labini et al 
al. |l99g|) and Z^ism 



1.6 — 2.3 (Elmegreen 



19981, Joyce et 
Falgarone 



first hints that the power-law relation between velocity 
dispersion and size is also valid for stellar populations 
and that it extends beyond the size of Giant Molecular 
Clouds. Several observations confirm that the hierarchical 
structure of kinematically cold media is not only present 
in Milky Way molecular clouds, but is also found in other 
systems and o n larger scales. For examples, Vogelaar & 
Wakker (1994) found perimeter- area correlations in high- 
velocity clouds; power-law power spectra of HI emission 
were found in the Sma ll and the L arge Magellanic Cloud 
by Sta nimirovic et al. ( 199£ , 2001 ), and Elmegreen et al. 
( 2000 ), respectively; measurements of the HI distribution 
in galaxies of the M81 cluster reveal fractal structures on 
the galaxy disc scale (Westpfahl et al. 19991 ). 

On cosmic scales, up to about 100 Mpc, matter is also 
hierarchically organized. A common feature of the ISM 



1996| , Combes |199^ ), respectively. Thus the precise value 
of the fractal dimension does not seem to be universal, but 
a range between 1 and 2 appears frequent. 

All this may suggest that a general scale free factor is 
mainly responsible for the matter distribution and the dy- 
namics of cosmic structures, galactic disks and molecular 
clouds. Only one factor appears to be dominant over all 
these scales, namely gravity. 

Gravo-thermal experiments on isolated systems show 
that typically two possible states are reached asymptoti- 
cally, a high energy homogeneous state and a low energy 
collap sed st ate, with a halo-core struc ture (Lynden-Bell & 
Wood |1968|, Hertel & Thirring |l97l|, Aronson & Hansen 



19721) 



Thus to produce more inhomogeneous structures self- 
gravitating systems must be open, such as be subjected to 
time dependent boundary conditions. On cosmic scales the 
Hubble flow represents a time dependent boundary con- 
dition, and develops lumpy structures. On galactic scales 
down to molecular cloud scales an energy flow, maintain- 
ing the system out of equilibrium, may be sustained by 
the shear-flow and small scale dissipation. Indeed, gravita- 
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tional instabilities convert directed kinetic energy (shear- 
flow) into thermal and turbulent motion (von Weizsacker 



1951, Goldreich & Lynden-Bell 1965). Turbulent motion 



may then transport the energy through the scales until it 
is dissipated away by radiation in molecular collisions and 
shocks. 

The lumpy distribution of matter reported by Toomre 
( |1990D and Toomre & Kalnajs (1991, hereafter TK) in 
local shearing-sheet experiments of disks reminds us of 
the ubiquitous inhomogeneous state of the ISM as well 
as the flocculent structures of many spirals. The rele- 
vance of these experiments for galaxies is supported by 
the recurring spirals found in slightly dissipative complete 
s elf-gr avitating disk simulations by Sellwood & Carlberg 
( 1985 ) and many others (e.g., Miller, Prendergast & Quirk 



197C ) . The TK models confirm that purely self-gravitating 



systems with time-dependent boundary conditions can 
produce very chaotic inhomogeneous structures. 

Here, in order to investigate in more detail the mat- 
ter distribution produced by self-gravitation, shear and 
dissipation we perform further such local shearing-sheet 
experiments. To check if the resulting structures reveal 
power-law relations we calculate the power-law indices of 
the mass-size and the velocity dispersion-size relation. To 
compare with earlier models, in particular with those of 
TK, we start with 2D simulations and extend then the 
model to 3 dimensions. This extension is important be- 
cause as soon as dense clumps develop in a disk with hor- 
izontal sizes comparable to or smaller than the supposed 
thickness of the disk, motion transverse to the plane must 
be strongly coupled to the motion in the plane, and the 
2D approximation is no longer valid. 

To obtain instructive models it is important not to in- 
clude too many ingredients. We are primarily concerned 
not with complex physical objects such as molecular 
clouds, but with processes. So our approach is not to in- 
clude a maximum of physical ingredients, but just the ones 
that appear as the most relevant. We want to check if self- 
gravitation in combination with time-dependent boundary 
conditions and a slight dissipation can produce and main- 
tain an inhomogeneous, lumpy and eventually self-similar 
structures, resembling those observed in galactic disks and 
molecular clouds. 

The considered scales are of the order of 0(1 — 10) kpc. 
Thus the transition regime between the molecular cloud 
scale (« 0.05 kpc) and the galactic disk scale (« 10 kpc) 
can be investigated. However, since the model is scale- 
free, we can draw conclusions beyond the here considered 
scales and thus eventually contribute to illustrate the scal- 
ing laws observed in sub- or extra-galactic structures. 

Preliminary results of our numerical experiments were 



presented in Huber & Pfenniger (1999, 2001). Since then 
we continued to improve our model and to collect more 
experience, which led to new insights with respect to the 
clustering simulation and scaling laws. In this paper we 
discuss in detail the model and the results. Similar studies 



In the next section, we justify the use of dissipative 
particles in order to model the dynamics of self-gravitating 
gas. The numerical model in presented in Sect. 3 and a 
pseudo-code is given in Appendix A. In Sect. 4 we discuss 
the methods used to analyze the structures resulting from 
our shearing box simulations. The results of the 2D and 
3D simulations are presented in Sect. 5. Finally, Sect. 6 is 
dedicated to discussing limitations of the models. 

2. Physical Gas Model 

2.1. Hierarchical Systems 

For a hierarchical self-similar structure, one can define a 



fragmentation efficiency (Scalo 1985) 
/ = ■qmL-i/mL , 



(1) 



have been presented by Semelin (Semelin 1999, Semelin & 
Combes pOOOl). 



where m^-i and are the mean masses of a fragment 
at level L—\ and L, respectively. The factor r] is the num- 
ber of fragments formed at each level. The fragmentation 
efficiency / indicates how much mass in a clump is con- 
centrated in subclumps. If / is not very high (< 95%), 
the smallest fragment masses become negligible after a 
few levels and a hierarchical description is less relevant. 
However, if / is very high, several iteration steps can be 
carried out and the bulk mass is still concentrated in the 
smallest subclumps. As long as the bulk mass is concen- 
trated in subclumps the interclump mass can be neglected. 
For convenience we call the smallest clumps for which the 
interclump medium can be neglected basic-clumps. If the 
level of the basic clumps is zero, a clump at level L is 
formed by 77^ basic clumps. 

Observations of the interstellar medium reveal a highly 
inhomogeneous and clumpy structure (see, eg.. Dame et 
al. |2000| , Tauber et al. |1991[ Storzer et al. pOOOj ). Moreover 
the structure is for a certain scale range hierarchical. 
Assuming that the size of particles is larger or equal to 
the size of the basic gas clumps, the structure of the in- 
terstellar medium can be described correctly down to the 
scale of the particles by the distribution of these parti- 
cles. A particle represents then the lowest resolvable level, 
while clumps at higher levels, i.e., at larger scales are rep- 
resented by an ensemble of particles. 

2.2. Dissipative Particles 

At the here considered scales, larger than tens of pc, the 
description of the dissipative processes taking place in the 
ISM is very complex and far from respecting the hypothe- 
ses allowing the full application of Navier-Stokes equa- 
tions. Thus the use of a traditional hydrodynamic code is 
in no way "better" than the simpler approach adopted by 
TK, where a simple small drag parameter is all what is 
introduced as dissipation. 

Indeed, we recall the following considerations: 

1.) Being long range, gravity breaks the fundamental as- 
sumption made in classical thermodynamics that interac- 
tions are short ranged. In turn, when gravity is sufficiently 
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strong (i.e., the Jeans' instability threshold is reached), 
supersonic chaotic motion is expected, as also systemat- 
ically reported for the interstellar medium. This means 
that no local pressure equilibrium is reached at the scales 
over which turbulence exists. Down to the smallest scale at 
which supersonic turbulence exists no local thermodynam- 
ical equilibrium can be established, and thus no equation 
of state can be defined. A basic assumption allowing to de- 
rive the usual Navier-Stokes equations of fluid dynamics is 
missing. Besides, numerous observational evidences indi- 
cate non-thermal cloud clumps. For instance, Beuther et 
al. (|200C| ) carried out multi-wavelength observations and 
compared the line ratios with radiation transport mod- 
els. They found that models based on the assumption of a 
local thermodynamical equilibrium (LTE) can not repro- 
duce the observed data set. Due to the lack of a LTE down 
to smallest scales, thermal physics appears as an inappro- 
priate tool to represent the statistical state of interstellar 
gas. 

2.) Fundamentally the Navier-Stokes fluid equations de- 
scribe a) the local conservation of mass and momentum, 
with b) additional constraints such as local smoothness 
of the quantities subject to differentiation, c) an equation 
of state for closing the moment equations derived from 
Boltzmann's equation, and d) phenomenological laws de- 
scribing viscous forces. While the mass and momentum 
conservation laws are likely to be adapted even for such 
a clumpy medium as the ISM, the other constraints do 
not. In this context the energy equation is little relevant 
(is not a constraint) if no control can be performed on ra- 
diative processes, which operate on very short time-scale 
in the cold ISM. Since large but clumpy entities such as 
molecular clouds have a mean-free path much larger than 
their size, the dynamics of such systems may as well, in 
the present state of understanding, be de scribe d by semi- 
collis ional, dissipative pa rticles (Brahic 1977 Pfenniger 
1998|) . Casoh & Combes ( |l982| ) studied the formation of 



giant molecular clouds through cloud collisions and coa- 
lescence in the molecular ring. They found that the en- 
semble of clouds never reaches a steady state. Thus they 
concluded that clouds are better described by particles 
than by a fluid. Another hint that the usual fluid equa- 
tions are not better adapted to describe the ISM than 
sticky particles is that the rings in barred galaxies are 
never reproduced by the former but easily by the latter 
(e.g., Schwarz|l984|). 



3.) In the ISM not only the cooling and heating processes 
are rapid with respect to global dynamics, but also the 
energy reservoir of global dynamics is much larger than 
the other available energy reservoirs represented by gas 
pressure, stellar radiation, cosmic rays, or magnetic fields 
(Pfenniger 1996). The virial theorem, expressing a bal- 
ance of negative and positive energy reservoirs, is a useful 
tool to order the importance of respective physical factors 
according to their quantitative values. Since the energy 
budget at the galactic scale is dominated by dynamics, 
to first order the system is well described by conserva- 



tive dynamics, and dissipative effects are of second order. 
In weakly dissipative systems the stable periodic orbits 
and fixed points of the conservative case are transformed 
into attractors, and chaotic orbits typically converge to- 
ward strange attractors with similar chaotic properties. 
Therefore one can naturally infer that the exact dissipa- 
tive force is irrelevant as long as it remains weak, since the 
long term behavior is an attractor. The dissipative pertur- 
bation is weak when during the time-scale of interest the 
energy dissipated is small with respect to the total energy 
of the system. Therefore in this regime it is not necessary 
to know precisely how energy is dissipated, any weak fac- 
tor le ads to the same attractors (see Pfenniger & Norman 



1990| for an extended discussion on the topic) . 



These considerations show that weakly dissipative par- 
ticles are a permissible method to study the dynamics of 
interstellar gas at sufficiently large scales. The mass and 
momentum conservation is granted by the equations of 
motion, and the weakly dissipative regime by a simple 
linear friction law. 

3. Numerical Model 

In previous studies of shearing sheet disks, the forces of the 
self-gravitating particles were computed by direct sum- 
mation. Instead, we show that by using a time-dependent 
affine coordinate system we can represent the shear-flow 
in periodic coordinates. Consequently we can increase 
the computation performance by calculating the self- 
gravitational forces with the popular FFT-convolution. 

3.1. Principle 

Here we explain the principle of the local model for the 
3D case. Ignoring all expressions with a z, yields the 2D 
case. 

In a local model of a disk, everything inside a box of 
a given size is simulated, and more distant regions in the 
plane are represented by replicas of the local box (Toomre 
& Kalnajs [l99l| . Wisdom & Tremaine |l988i Salo [19951) . 
The global galactic disk attraction made by components 
such as stars or dark halo not included in the local box is 
1.) cancelled to zeroth order by adopting a rotating frame, 
and 2.) corrected to first order by the linear terms in the 
epicycle k and vertical frequency v (Binney & Tremaine 
1991 ). 

In the same spirit as TK, the matter in the box has an 
undefined mass composition with a slight dissipation. For 
a normal spiral each particle may be considered as a mix- 
ture of stellar mass and gas, with mean weak dissipation. 

The origin of the particle coordinates in the local box is 
a reference point that moves on a circular orbit at distance 
TZq from the galactic center with the orbital frequency 
r^o = ^{T^o)- In their model, TK used a rotating Cartesian 
coordinate system. The horizontal particle positions are 
then given hy x — TZ — TZq, y = TZo{0 — rioO ^nd the 
vertical location by z. If x^y,z <^ TZq, the orbital motion 
of the particles is determined by Hill's approximation of 
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Newton's equations of motion (Hill 1878). In the present 
context, they read: 



y + 2Qox = Fy 

z = -v'^z + F, 



(2) 



where Aq = ~ ^TZo{di} / dTVj-jZo ^^e Oort constant of dif- 
ferential rotation. F^ , Fy and Fz are local forces due to 
the self-gravitating particles, that should be small with 
respect to the global force field. Like TK and Griv et al. 
( |1999D we use these equations also for simulation zones, 
where x, y, 2: <C TZq is not valid for the most part, i.e., for 
galactic disc scales, meaning that non-linear higher order 
effects are not taken into account. However, since much 
of the gravitational force in any wavy disturbance stems 
from the nearest particles (TK, Julian & Toomre 1965, 
Toomre 1964), the conclusions of the model should be rel- 
evant for galactic disks, despite the violation of the lin- 
earity hypothesis. Indeed, the swing amplification theory, 
whose applicability to spirals has been well established, is 
based on the same assumptions. 

In a Cartesian coordinate system (x, y, z) the posi- 
tions of the rectangular boxes (local simulation box and 
replicas) change with time due to the differential rota- 
tion dflo/dx. The differential rotation causes a shear-flow, 
which reads for a flat rotation curve, y = flox. Thus a 
particle at {x, y, z) has images at (x + nLx, y — nLxilot + 
mLy, z), where n and m are integers (Wisdom & Tremaine 
1988). Lx, Ly and Lz are the sides of the local box. An 
initially (t = 0) periodic arrangement of the boxes rel- 
ative to a fixed Cartesian coordinate system can not be 
maintained (see Salo 1995 ). As a consequence the forces 
of the self-gravitating particles have been determined in 
previous simulations by direct summation with upper and 
lower cut-offs, meaning that the computation time for N 
particles is proportional to N"^. 

However we can improve the performance by com- 
puting the forces in the Fourier space with the convolu- 
tion method and the FFT algorithm (Press et al. 1986|) . 
Thereby the potential computation time is reduced to be 
proportional to Nclog{Nc), where Nc is the number of 
cells, taken here as proportional to the number of particles. 
The FFT approach requires a system spatially isolated 
and/or periodic at all times. Here the system, represent- 
ing the local dynamics of a disk, is isolated in z-direction. 
In the X — y— plane the system is periodic, but only on 
affine coordinate systems whose pitch angles change peri- 
odically. 

The dark box in Fig. la represents the local box in a 
Cartesian coordinate system (solid mesh). In the initial 
state a certain local particle (star in the dark box) and 
its replicas (stars outside the box) are periodic relative to 
the Cartesian coordinate system. But then the particles 
are shifted by the shear and the periodicity relative to 
a rectangular coordinate system is lost. However because 
the shear is linear in x there is for all times an affine coor- 
dinate system {x',y',z') on which the system is periodic. 



Thus we modify our initially rectangular coordinate sys- 
tem with a time dependent pitch angle. The solid mesh in 
Fig. la-c represents an afhne coordinate system in which 
the periodicity of the system is maintained. Its pitch angle 
is, ab < 0, for all times. Thus we call this coordinate sys- 
tem for convenience the backward mesh. The inclination 
of the backward mesh ai, is determine by the shear. 



at — tana;, = [{—ilcit) inod{Ly / 



(3) 



Fig. Ic shows the system at t = Ly/ {Lx^q). We can see 
that the periodic arrangement of the particle images cor- 
responds to those in Fig. la. Consequently the system is 
again periodic on a rectangular coordinate system and we 
can replace the backward mesh in Fig. Ic with the one 
in Fig. la. Thus the inclination of the afhne coordinate 
system jumps at i = Ly/(Lx^o) from at = Ly/Lx to 
ab = 0. This accounts for the modulo function in Eq. (||). 
Thus, if Lx/{LyD,o) is a multiple of the time-step, only 
a finite number of affine coordinate systems is necessary. 
As a consequence the corresponding kernels must be com- 
puted only once at the beginning of the simulation and 
stored for subsequent use. 

For the computation of the forces of the self-gravitating 
particles one coordinate system in which the matter dis- 
tribution is periodic at all times would in principle be 
enough (e.g., the backward mesh). However, in order to 
avoid discontinuities in the force field when the inclina- 
tion of the coordinate system ab jumps back to zero, we 
compute the forces additionally in a second affine coordi- 
nate system, in which the system is periodic as well. The 
dashed mesh in Fig. la-c represents this second coordinate 
system. Because its pitch angle is always, a/ > 0, we call 
it the forward mesh. The inclination of the forward mesh 
a/ can be deduced from those of the backward mesh by: 



2j: = tana/ ^ at + (Ly/Lx) , 



(4) 



The light box in Fig. p.l| is the local computation box of 
the forward mesh. After the computation of the forces F' 
in both coordinate systems, we add them with weighting 
factors in order to soften the effects of the abrupt tran- 
sition a.t t ~ Ly/{Lx^o) on the force field. The forces 
computed in the forward and the backward mesh are 
and F'j:, respectively. Before adding the forces with the 
corresponding weighting factors, we transform them to 



Cartesian coordinates, F^ 



Fb,F'j: Ff. The single 



components of the forces are transformed as follows: 



F — 

F — 
F^,z = 



F' 

F' -+ 
F' 



aF! 



(5) 



where i = b, f for the backward, respectively for the 
forward mesh. Then the forces are weighted and added, 
F = bFt, + fFf. The weighting factors b and / are nor- 
malized (6 + / = 1) and proportional to the mesh incli- 
nation, b = —aiiLx/ Ly. Because forces are additive such 
a weighted force summation is permissible. The forces F 
correspond now to those in Eq. (S). That is, the inclined 
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forward mesh Orbital Motion 




Fig. la-c. The dotted frame represents a section of the disk, being infinite in two directions, seen from above. The two meshes are 
affine coordinate systems on which the mass distribution is periodic. The dark and the hght box represent the local computation 
box in affine coordinates, i.e., in the forward and the backward mesh respectively, a: The initial state of the two meshes {t = 0). 
The inclinations of the meshes are ai, — and a/ = (Ly/Lx)- Thus the corresponding weighting factors are b = 1 and / = 0, 
respectively. Consequently, the forces in the Cartesian coordinates are for this situation, F = Ft. Below the two meshes, the 
pitch angles of the affine coordinate systems are indicated, q/ and at, respectively. The angle between the two meshes remains 
the same at all times (a/ + aj, = a,Tcta.n{Ly / Lx))- b: The meshes and the weighting factors ai t — Ly/{2Lx^lo)- It is valid, 
Ub = —Ly/{2Lx) = —a/ and thus F = Ff,/2 + Ff /2. c: When the meshes reach these inclinations (t = Ly / [LxQ,o)), they jump 
back to the positions shown in (a) and the process starts again without introducing discontinuities in the dynamics. 



coordinate system are only used to compute the forces of 
the self- gravitating particles with the convolution method; 
then they are transformed to a Cartesian coordinate sys- 
tem. The evolution of the system in the Cartesian coordi- 
nate system is given by Eq. (||). 

The weighting described above, not only softens the 
effect of the abrupt change in time of the pitch angles, 
but also minimizes asymmetry effects due to the mesh in- 
clinations. Asymmetry effects disappear for example com- 



pletely when the inclination of one of the meshes is zero 
or when both meshes have the same inclination. In the 
first case (Fig. la) the weighting factor of the uninclined 
mesh is one and thus the forces are computed exclusively 
in the rectangular coordinate system. In the second case 
(Fig. lb) the asymmetry effects in both inclined coordi- 
nate systems cancel each other out. 

In an inclined coordinate system the gradient V de- 
pends on the inclination. This must be taken into account 
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by the calculation of and F'^ . The Euler-Lagrange equa- 
tions yield then the forces in an afRne coordinate systems, 



, _ (9$ a$ 



F' . 



where i = {6, /} for the backward, respectively for the 
forward mesh. 



3.2. Canonical Equations 



Pfenniger & Friedli (|1993| ) shown that the use of a leap- 
frog finite difference approximation of Newton's equations 
in a rotating reference frame with non-canonical variables 
lead to instability ("complex instability") in the sense of 
von Neumann. This is not the case when canonical vari- 
ables are used, then the stability or instability character is 
conserved between the leap-frog algorithm and the orbits. 
Therefore our model uses these equations. The canonical 
equations of motion with the momenta {px,Py,Pz} are. 



(7) 



z — 
and 

Px = 
Py - 
Pz = 



Pz 



nl)x 



Fx + 
F, - 



iloPy 

^OPx 



(8) 



These equations are invariant under the linear trans- 
formation 



X 

y 

Px 
Py 



X 

y 

Px 
Py 



2Ao t kx 
2Af)Q,Qtkx 
(r^o - 2Ao)kx 



^Oky 



(9) 



where kx and ky are arbitrary numbers. Thus, whenever 
a particle leaves the local box Lx x Ly x in the x or y- 
direction and its image enters somewhere on the opposite 
side (in the affine meshes the image enters exactly at the 
opposite face), we also have to transform the canonical 
momenta and their time derivatives correspondingly to 
the rules given above. 

3.3. Kernel 

For the 2D simulations we use an isotropic interaction po- 
tential. However, in order to resolve the flat disk vertically 
an anisotropic kernel is necessary due to computational 
limits. Thus most of the 3D simulations are carried out 
with an anisotropic kernel having the form of a paral- 
lelepiped. 



3.3.1. Isotropic Kernel 

In affine coordinates the softened isotropic interaction po- 
tential has the form. 



(6) $ 



( 1 fn {l+a^)x^+y^+2axy 
27 ' ^ ^ 



r < e 



r > e , 



(10) 



where a is the mesh inclination and e is the softening 
length. The advantage over a Plummer potential, used 
by TK and many others, is that this potential become a 
correct l/r gravitational potential beyond the softening 
length. Thus there is no sum up of small errors of the 
gravitational force due to the many dista nt part icles as in 
the case of a Plummer potential (Dehnen 2000 ) . 



3.3.2. Anisotropic Kernel 

The simulation box, representing local dynamics of a disk 
galaxy on the kpc scale, is rather flat {L^ <^ Lx,Ly). Thus 
our 3D-model needs an anisotropic force resolution and 
consequently an anisotropic kernel. This will be explained 
more exactly in the following. To calculate the forces of the 
self-gravitating particles we use a particle-mesh method. 
This method consists of three steps. First, the particle 
masses are assigned to the nodes of a mesh, which we 
call simulation mesh. We do this in accordance with the 
cloud-in-cell (CIC) scheme (see, e.g., Hockney & Eastwood 
1981). The masses at the nodes of the simulation mesh 
can be considered as new particles representing the mass 
distribution of the original particles. Second, the forces for 
the new particles are calculated on the simulation mesh 
node s via the convolution method (Hockney & Eastwood 
198l|) : 



^ijk — KijkPijk 



$ = $ 



ijk ' 



(11) 



(12) 



where ~ and ~ are the Fourier and the inverse Fourier 
transform, respectively. Kijk is the kernel and pijk is 
the mass density at a simulation mesh node r'--i.. If the 



mesh has Nx x Ny x N. 



ijk' 

,Nx 



nodes then i — 1, . . . , iv^,, ; j 
1, . . . ,Ny] k — 1, . . . ,Nz. The apostrophe ' indicates that 
the mesh is defined in inclined coordinates. 

The forces can now be calculated via Eqs. d). Finally 
the forces are interpolated at the original particle posi- 
tions. 

In order to avoid singularities and to approximate bet- 
ter physical objects self-gravitating particle are consid- 
ered to have an extent and consequently a softened poten- 



tial. Pfenniger & Friedli (1993) optimized the softening by 
adapting the particle extent as well as possible to the cell 
shape of the simulation mesh. In their polar-mesh simu- 
lations they approximated the cell shapes with a uniform 
triaxial ellipsoid. Here we adopt the particle extent to the 
cell shape as well. At a given time the cell shapes are all 
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simulation mesh 
new particle mesh 



a: pitch angle 



-J 



a 



Ix 



Ix 




ly 



Ly 



Lx 

Fig. 2. 2D representation of the simulation mesh and the 
new particle mesh depicting the discrete particle realiza- 
tion. The star indicates the origin. To calculate the kernel 
Kijk at the node (j, j, k) of the simulation mesh, one has 
to sum up over all mass points. The mass points are rep- 
resented by dots in the cell centers of the new particle 
mesh. 

identical, typically, because of the shear, a non-rectangular 
parallelepiped. In the orthogonal case the corresponding 
analytical form of the potential is known (McMillan 1958| ) . 
The analytical expression of such a potential is however 
quite cumbersome. Moreover we need also to describe the 
non-rectangular case. Thus we use a discrete realization 
of the particle mass. To this end we distribute the mass 
of the new particles over a refined discrete mesh having 
the same size as a cell of the simulation mesh. We call 
the refined mesh, new particle mesh. Simulation mesh and 
new particle mesh are shown in Fig. |^. To calculate the 
kernel Kijk at position r^^j, one has to sum up over all 
mass points of the discrete particle realization. 



ijk 



7V„ N^, Af„ 

EEE 

u—l v—1 w—l 



1 



ijk uvw 



(13) 



where x Ny x N^^ is the number of mass points rep- 
resenting a particle and r^„^ are the cell-centers of the 



new particle mesh (see Fig. g). The positions are given 
in affinc coordinates r' = {x',y',z'). In order to calculate 
the kernel the following coordinate transformation is thus 
carried out, 



y ax , 



(14) 



where a is the the inclination of the afhne coordinates. 
The inclination is fixed by the pitch angle, a = tana. 
Consequently the denominator in Eq. (O) has the form 



ijk uvw I 

((1 + a')ix', - x'J' + {yr - y'^f + (z^ - z'^ + (15) 
2a{x\-x'M,-v'.)f 



Nl/2 



Since the kernel needs to be evaluated only once for ev- 
ery possible inclination the cost of this procedure remains 
negligible. 

Because the cell size of the simulation mesh determines 
the particle extent, the softening is automatically fixed by 
the choice of the simulation mesh. 

It is important that the origin of the kernel represents 
the center of the simulation box in order to avoid a non- 
zero temporal mean velocity of the center of mass, which is 
introduced by an asymmetric description of the centrifugal 
force. Thus the positions r'^^^ = (x\,y'pz'^ are fixed as 
follows: 



(16) 



where lx x ly x Iz is the size of a simulation mesh cell. 
The positions representing the discrete mass distribu- 
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2 / ^ 


U = 


1, 








V = 


1, 


■ ■,N, 




^ N^ + l)l 
2 


W = 


1, 





(17) 



where l^xlyX is the cell size of the new particle mesh. 

The system is not periodic in the z-direction. To sup- 
press the images introduced by the FFT we use th e clas - 
sical doubling-up procedure (Hockney & Eastwood 1981 ) , 
which by doubling the size of the mesh over which the 
FFT must be performed exactly cancels all the images. 
Thus only the lower half of the entire mesh is relevant, 
i.e., only particles inside — Lz/2 < z < are active and 
particles leaving this zone are considered as escaped. 

3.4. Friction 

In the ISM the collisional rate must depend on the clump- 
ing state, which must depend on the dissipation rate. 
Consequently, we expect a complex dependence of drag 
coefficients and mass density. However, as explained in 
Sect. 2, at the kpc scale the physics is dominated by the 
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conservative gravitational dynamics and its concrete be- 
havior should be weakly dependent on the particular dis- 
sipative factors, since dissipation mainly acts to ensure 
the convergence of the system toward the attractors de- 
termined by the conservative part of the system. Thus, 
and following TK, as dissipation factors we adopt linear 
friction terms, which should be weak in order to remain 
quasi-Hamiltonian (Pfenniger & Norman 199C| ). 

Yet the collisional properties of the interstellar medium 
can be expected to differ along or transverse to the plane. 
To minimize the number of free parameters, we retain only 
two friction coefficients. The linear friction terms —C^i 
and —CzZ added to the radial respectively to the vertical 
forces {Fx,Fz) in Eq. (||). There is no azimuthal friction 
in order to be consistent with a global angular momentum 
conservation. 



3.5. Scaling, Units, Parameters 

In order to fix a scale, the origin of our local model is 
located at a distance of TZq = 8 kpc from the galactic cen- 
ter and rotates with an orbital frequency of fio = 
where 6*0 = 210 km/s. We assume for the general case a 
flat rotation curve. Moreover, we assume that the active 
disk has a surface density of Sq = 100 Mq/pc^. 

As usual in local shear models of galactic disks the 
linear measure is the critical wavelength, i.e., the longest 
unstable wavelength in a zero-pressure disk, 



A. 



(18) 



The critical wavelength defines the scale for which the 
theory of swing amplification predicts the strongest re- 
sponse (Toomre 1981 , Julian & Toomre 1965| ). For a flat 
rotation curve the epicyclic frequency is, k = \/2f^o and 
consequently the critical wavelength scales, with the pa- 
rameter values indicated above, to 1 Adit = 12.32 kpc. The 
disk scale height zq is then 0.024 Adit- However, the equa- 
tions of motion are scale free and the model can, with an 
appropriate choice of the parameters, be rescaled at will. 

The friction coefficients Cx and Cz of the friction terms 
—Cxi and —CzZ are in this work indicated in units of 
1 / Tosc , where Tqsc is the period of the unforced epicyclic 
motion. The cooling times of the radial and the verti- 
cal damping are thus tcooi.s = 1/Cx Tosc and tcooLz ^ 
1/Cz Tosc- For all models presented here, Tosc < ^cooi ap- 
phes. 

The time-step has to meet the following conditions: 



At < 0.1 min{Zi/fTi} , i — {x,y,z} 



At = - 



1 Lx 



(19) 
(20) 



where ai is the initial velocity dispersion ellipsoid, k is the 
cell size and k is an integer. According to Eq. (^ the evo- 
lution of the inclination of the backward grid is periodic 
with period T = Lx/{LyQo). The second condition guar- 
antees that this period is a multiple of the time-step. Thus 



the number of possible grid inclinations and consequently 
the number of kernels is finite. In order to satisfy the above 
conditions the time-step is computed in two steps: 



k 
At 



lOLx 



Lyilo min{/i/(Tj 
Lx 



(21) 
(22) 



where [] means to round to the next higher integer. 

TK calculated the forces of the self-gravitating parti- 
cles with direct summation. Thus they had to introduce 
an upper cutoff in order to limit the computational ex- 
penditure, meaning that beyond a certain separation the 
particles lost their mutual gravitational interaction. Their 
separation cutoff was equal to four times the softening 
length. This limited the dynamical range of gravity to 0.6 
dex. They argued that a cutoff at larger separations did 
not affect the resulting structures. Thanks to the higher 
performance of the convolution method we can extend 
the dynamical range without increasing the computation 
time. This may be important in view of a self-similar mat- 
ter organization in self-gravitating systems. 

The parameters characterizing the model of TK are 
indicated in Table |^ and ^. They carried out numerical 
shearing-sheet experiments for different particle densities 
n. Their friction coefRcient is a function of the particle 
density Cx = (3.5 x 10~^)/n t~J.. In order to extend this 
study and to explore the resulting structures in depen- 
dence of the different parameters, we realize different ver- 
sions of the shearing box model. These model versions are 
characterized by different parameter sets which fix the size 
of the simulation zone, the resolution, the particle density 
etc.. For convenience we call these model versions in the 
following models. That is, a model denotes in the following 
a version of the shearing box model which is determined 
by a specific parameter set. The parameters of the models 
are indicated in Table || and ^. 

To be able to do some statistics of structures produced 
on scales with strongest swing amplification response we 
perform, like TK, simulations with a quite large simula- 
tion zone. Model 1-11 have such a large simulation box 
resp. simulation sheet in the 2D case, Lx x Ly{xLz) = 
6 X 6(x0.8) Aj?,it. However since the dynamical range is 
limited due to computational limits we perform also sim- 
ulations for a smaller local box, in order to resolve smaller 
scales. Therefore the simulation box of model 12-19 are 
reduced to 1.8 x 1.8 x 0.8 

The time-step depends on the mesh resolution. The 
mesh resolution is fixed by the number of particles and 
the size of the simulation zone. The computation time of 
a simulation depends thus on the particle density n. We 
increase n with respect to previous models based on direct 
force calculation up to a factor 30 and are furthermore 
able to perform the simulations in 3D. The code has been 
written in Matlab for its ease of use, but clearly a compiled 
language program would greatly improve its speed and 
memory usage. A pseudo code is given in Appendix A. 
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Model 


Lx X Ly 


n 


Dynamical range 


Potential/ 


# Dimensions 




[^crit] 




[dex] 


Softening 




TK 


6.0 X 8.0 


100-1200 


0.6 


Plummer 


2 



Table 1. Parameters characterizing the model of TK. Indicated are the size of the simulation zone, the number 
density of particles (surface density), the dynamical range, the gravitational potential of the particles and the number 
of dimensions. 



Model 


a/10-'^ [I/Tosc] 


iV 


£ [Acrit] 


k/SIo 




TK 


3.5/n 


4800 - 57000 


0.20 


1.4 


0.5 



Table 2. Parameters of the TK model. The friction coefficient Cx is a function of the particle density n. N is the 
particle number, e is the softening length of the Plummer potential. The epicycle frequency k and Oort's constant 
are indicated in units of fJo- 



Model 


Lx X ILy X Ijz 




n 


Dynam. range [dex] 


Potential/ 


# Dim. 


Var 




[Acrit] 


[Acrit] 




plane 


vertically 


Softening 








1 


6.0 X 6.0 A^Ht 


0.023 X 0.023 A^Ht 


1820 


1.5 




Isotropic 


2 


Cx 




2 


6.0 X 6.0 Acrit 


0.023 X 0.023 AcHt 


1820 


1.5-2.5 




Isotropic 


2 


e 




3 


6.0 X 6.0 X 0.8 


0.188 X 0.188 X 0.013 


910 


1.3 


0.5 


Isotropic 


3 






4 


6.0 X 6.0 X 0.8 


0.094 X 0.094 X 0.013 


3640 


Var 


Var 


Isotropic 


3 


e 




5 


6.0 X 6.0 X 0.8 


0.188 X 0.188 X 0.013 


910 


1.5 


1.8 


Anisotropic 


3 


Cx , 


V = 0.3 


6 


6.0 X 6.0 X 0.8 


0.188 X 0.188 X 0.013 


910 


1.5 


1.8 


Anisotropic 


3 


Cx , 


V = 3.0 


7 


6.0 X 6.0 X 0.8 


0.188 X 0.188 X 0.013 


910 


1.5 


1.8 


Anisotropic 


3 


a 




8 


6.0 X 6.0 X 0.8 


0.188 X 0.188 X 0.013 


910 


1.5 


1.8 


Anisotropic 


3 


V 




9 


6.0 X 6.0 X 0.8 


0.094 X 0.094 X 0.013 


3640 


1.8 


1.8 


Anisotropic 


3 


Cx 1 


V = 0.3 


10 


6.0 X 6.0 X 0.8 


0.094 X 0.094 X 0.013 


3640 


1.8 


1.8 


Anisotropic 


3 


Cx 1 


V = 3.0 


11 


6.0 X 6.0 X 0.8 


0.094 X 0.094 X 0.013 


Var 


1.8 


1.8 


Anisotropic 


3 


N 




12 


1.8 X 1.8 X 0.8 


0.056 X 0.056 X 0.013 


10100 


1.5 


1.8 


Anisotropic 


3 


Cx 1 


V = 0.3 


13 


1.8 X 1.8 X 0.8 


0.056 X 0.056 X 0.013 


10100 


1.5 


1.8 


Anisotropic 


3 


Cx , 


V = 3.0 


14 


1.8 X 1.8 X 0.8 


0.056 X 0.056 X 0.013 


10100 


1.5 


1.8 


Anisotropic 


3 






15 


1.8 X 1.8 X 0.8 


0.056 X 0.056 X 0.013 


10100 


1.5 


1.8 


Anisotropic 


3 


u 




16 


1.8 X 1.8 X 0.8 


0.028 X 0.028 X 0.013 


40450 


1.8 


1.8 


Anisotropic 


3 


Cx , 


V = 0.3 


17 


1.8 X 1.8 X 0.8 


0.028 X 0.028 X 0.013 


40450 


1.8 


1.8 


Anisotropic 


3 


Cx , 


V = 3.0 


18 


1.8 X 1.8 X 0.8 


0.028 X 0.028 X 0.013 


40450 


1.8 


1.8 


Anisotropic 


3 


a 




19 


1.8 X 1.8 X 0.8 


0.028 X 0.028 X 0.013 


40450 


1.8 


1.8 


Anisotropic 


3 


Ao 





Table 3. We use 18 models to explore the different parameters. Besides the parameters presented in Table |l| for the 
TK model, we indicate here the mesh resolution IxXlyX Iz and the dynamical range vertical to the plane. Var indicates 
the parameter, altered from run to run. The resulting structure are then explored as a function of this parameter. The 
gravitational particle potential is either isotropic or it is deduced from the discrete particle representation described 
in Sect. |3.3.2| 



3.6. Code Testing 

In order to check our code we carry out two-body simula- 
tions and compare the results with analytical solutions of 
the Kepler problem. We use a non-rotating inertial frame, 
thus ^Iq = Aq — — in the equations of motion (Eq. 
and (H)). However we compute the forces on the time 
dependent afhne coordinate systems. Thus flo — Oq/TZq 
in Eq. (H), where the inclination angles are determined. 
That is, we calculate the forces for a non-rotating iso- 
lated system with the help of "shearing Fourier meshes" . 
Furthermore we use the anisotropic kernel described in 
3.3.2| . 



Sect. 



we test our code for different resolutions I. The particle 
extension is fixed by the anisotropic kernel and is equal to 
the resolution. 

In the initial state the velocities of the two particles 
are chosen, in the way that they move on circular orbits. 
Accordingly to theory the following holds at each time: 



Ar 



r - ro = 



= 



(23) 



The vertical resolution of all simulations presented in 
this work is 1^ — 0.013 Adit, but the resolution in the 



plane I = = ly depends on the particle number. Thus T — 2n\ - 



where r is the relative particle distance at i > and rg is 
the initial particle distance at t = 0. rem is the center of 
mass. The orbital period for a particle with mass m is, 



(24) 




10 



D. Huber, D. Pfenniger: Lumpy Structures in Self-Gravitating Disks 



Model 


Ci/10 

[l/Tosc] 


Cz/10 

[l/ToscJ 


N 


e 

r \ 1 
[Acrit] 






Ao/ih 


Var 




1 


Var 


- 


65520 


0.2 


- 


1.4 


0.5 


Cx = 


40 - 210 X 10"^ 


2 


100 


- 


65520 


Var 


- 


1.4 


0.5 


e = 


0.02 - 0.3 


3 


Var 


0.7 


32760 


0.3 


0.3 


1.4 


0.5 


Cx = 


40 - 280 X 10"^ 


4 


100 


0.7 


131040 


Var 


0.3 


1.4 


0.5 


e = 


0.1 - 0.4 


5 


Var 


0.7 


32760 


- 


0.3 


1.4 


0.5 


Cx 


40 - 280 X 10'^ 


6 


Var 


0.7 


32760 


- 


3.0 


1.4 


0.5 


' = 


70 - 280 X 10"^ 


7 


140 


Var 


32760 


- 


0.3 


1.4 


0.5 


Cz = 


0.04 - 40 X 10"^ 


8 


140 


0.7 


32760 


- 


Var 


1.4 


0.5 


u = 


0.0-6.4 


9 


Var 


0.7 


131040 


- 


0.3 


1.4 


0.5 


Cx — 


40 - 210 X 10"^ 


10 


Var 


0.7 


131040 


- 


3.0 


1.4 


0.5 


Cx — 


- 120 X 10"^ 


11 


70 


0.7 


Var 


- 


0.3 


1.4 


0.5 


N = 


16000 - 128000 


12 


Var 


0.7 


32720 




0.3 


1.4 


0.5 


Cx — 


10 - 70 X 10"" 


13 


Var 


0.7 


32720 




3.0 


1.4 


0.5 


Cx — 


30 - 50 X 10"^ 


14 


50 


Var 


32720 




0.3 


1.4 


0.5 


c. = 


0.04 - 40 X 10"^ 


15 


50 


0.7 


32720 




Var 


1.4 


0.5 


1/ — 


0.0-6.0 


16 


Var 


0.7 


131050 




0.3 


1.4 


0.5 


Cx 


20 - 40 X 10"^ 


17 


Var 


0.7 


131050 




3.0 


1.4 


0.5 




20 - 40 X 10"^ 


18 


30 


Var 


131050 




0.3 


1.4 


0.5 


' = 


0.04 - 40 X 10"^ 


19 


30 


0.7 


131050 




0.3 


Var 


Var 


K = 


1.4, 1.7; = 0.5, 0.25 



Table 4. Parameters of model 1-19. Contrary to the TK model Cx 
of the extension to three dimensions, a vertical friction coefficient Cz 
for which we explore the models is indicated in the last column. 



is a free parameter. Moreover we have, because 
and a vertical frequency v. The parameter range 




0.02 



0.015 



-12.5 



0.01 



< 



0.005 



^ crit' 

Fig. 3. Relative errors, resulting from the simulation of 
two bodies on circular orbits, as a function of the resolu- 
tion. The simulations are performed for one period, which 
corresponds to « 4 galactic rotations. Crosses, left ordi- 
nate: The maximal relative error of the particle distance 
as a function of the resolution. The resolution is indicated 
in units of Acrit- Stars, right ordinate: The relative error 
of the orbital period. 



Particle mass, distance and velocities chosen for these tests 
yield a period of T « 4 galactic rotations. 

These theoretical results are compared with the exper- 
imental results, i.e., with those resulting from our simu- 
lations. The code errors, arising from this comparison are 
shown in Fig. ||and0 for different resolutions. The resolu- 
tion is indicated in units of Acrit- During an orbital period 
Ar/j'o oscillates around zero. The maximal error of the 



-13.5 



-14.5 




Fig. 4. The maximal acceleration of the center of mass 
during the simulation of two bodies on a circular orbit as a 
function of the resolution. The simulations are performed 
for one orbital period. 



particle trajectory computed with the numerical model 
are then equal to the amplitude of this oscillation. The 
amplitude (Ar/ro)max is plotted in Fig. ^. In this figure 
the relative error of the orbital period AT/T is indicated 
as well. Fig. ^ reveals the acceleration of the center of 
mass rem- The here presented errors must be considered 
as upper limits, because our particles are not point-like 
but have an final extension, which is, as mentioned above, 
equal the resolution I. In Table |^ we indicate the ratio l/r^ 
for the different resolutions. 
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I [Acrit] 


l/ro 


0.188 
0.094 
0.047 


0.094 
0.047 
0.024 



Table 5. A small contribution to the deviation from the 
theoretical trajectories is due to the extension of our test 
particles. In this Table we indicate the ratio of particle 
extension and separation I / rg for the different resolutions 



3.7. Initial Conditions 

Because we are interested in the secular time behavior of 
the galaxy disk, the simulations are performed for t = 10 
galactic rotations. In the initial state at i = the par- 
ticles are distributed uniformly in the x-y-plane. In the 
z-direction the particle distribution follows an isothermal 
law 



p cx sech^(z/zo) , 



(25) 



where p is the density and zq is the disk scale height. The 
velocities at t = are determined by the shear 



i = 

y = -2Aox 

i = 

and the Schwarzschild velocity ellipsoid 
3.36GSoQ 

C^x — 

K 

_ CTxK 



(26) 



(27) 



where th e Saf ronov-Too mre s tability criterion is Q > 1 
(Toomre 1964 , Safronov 1960 ). This velocity distribution 
is a permissible assumption, because we represent the gas 
by dissipative particles. 

4. Structure Analysis 

In order to characterize the structures resulting from the 
shearing box experiments, we determine the mass-size re- 
lation and the velocity dispersion-size relation. 

4.1. Mass-Size Relation 

We choose randomly a set of particles with distances 
r < L/4 from the center of the simulation box. The posi- 
tions of the particles in the set are restricted to r < L/4 
in order to avoid boundary effects in the analysis of the 
mass-size relation. We will refer to that at the end of this 
subsection. For each particle in the set we count the num- 
ber of neighboring particles N{R) inside a certain radius 
(all particles in the simulation zone are considered as pos- 



sible neighbors). If we repeat this for other values of R we 
can find the structure dimension D[R) via 



D[R) 



d\TL{N) 



(28) 



d\n{R) 

where R denotes the scale. The mass-size relation is then 

N{R) cx M{R) (X R^'-^^ . (29) 

If the structure dimension is independent of the scale, D = 
Df, i.e., if D is constant or oscillates around a mean value , 
then the mass-size relation is a power- law (Semelin ^000| ) 



M cx i?^ 



(30) 



If furthermo re Df is non-integer, the structure is fractal 
(Mandelbrot |1982D. We will check i{D = Df holds for the 



structures resulting from the shearing box experiments. 
However one has to take into account that these struc- 
tures can never correspond to an idealized mathematical 
set, generated by means of an infinite number of levels. 
Rather they are the result of a finite simulation, model- 
ing a finite physical system. Thus the structures resulting 
from our experiments can never be fractal beyond a lower 
and upper cutoff. An upper scale limit due to the numer- 
ical model is given by the size of the simulation box in 
the a;-y-plane. On this scale the system becomes periodic, 
meaning that it can not be fractal. To avoid boundary ef- 
fects as much as possible, we consider in our analysis only 
particles inside the simulation box. Therefore the frac- 
tal dimensions are only calculated for scales R < L/A. 
Moreover we determine the number of neighboring parti- 
cles only for particles with a distance r < L/A from the 
center of the simulation box. Consequently even for par- 
ticles at r = L/4 all neighboring particles are inside the 
simulation box. A lower scale limit is due to the finite res- 
olution of the simulation mesh. If the mesh cells have the 



size 



and I = Ix = ly > Iz then we do not expect 



to model correctly fractal structures below 21. 

If however the structure dimension depends on the 
scale D = D{R), the structure dimension may simply be 
regarded as a statistical measure describing the clumpy- 
ness on the corresponding scale. 

4.2. Velocity Dispersion-Size Relation 

There is observational evidence for Larson's law 



(31) 



on scales 0(0.1) - 0(100) pc with 0.3 <, Sl ^ 0.5 (e.g. 



Larson 1981, Scalo |1985| , Falgarone & Perauh |1987t Myers 
& Goodman 1983[ ). This power-law relation seems to ex- 
tend beyond the 100 pc scale (Larson 1978 ). 

In order to check if our model can reproduce these 
or similar velocity-size correlations on the kpc scale, we 
determine the velocity dispersion-size relations for the re- 
sulting structures. We use the same approach as for the 
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determination of the fractal dimension, but now we cal- 
culate the velocity dispersion a of the particles inside a 
certain radius R. Then 



5{R) 



diner 
dlni? 



(R) 



(32) 



What we said about the structure dimension applies also 
for 5. If 5 is constant or oscillates around a mean value over 
a certain scale range it may be regarded as the power-law 
index of Larson's law on the kpc scale, ^ = 5^. If however 
6 depends on the scale, 5 = S{R), it may be considered as 
a statistical measure determining the velocity correlation. 

5. Results 

5.1. 2D Simulations 

To compare with earlier simulations, we start with 2D 
shearing sheet simulations (model 1, 2 in Table 3, 4). 
The structures resulting from these simulations depend 
mainly on the relative strength of the competing gravita- 
tional and dissipation processes. Even without dissipation 
filamentary structures are already developed after w 1/2 
galactic rotation. Yet, gravitational instabilities lead to a 
conversion of bulk kinetic energy (shear-flow) into ran- 
dom thermal motion. In this way the disk is heated up. 
Thus, if there is no or a too weak dissipation the ini- 
tially arised filamentary structures are not maintained 
and smear out. However, with an appropriate dissipation 
strength, filamentary structures can be maintained in a 
statistical equilibrium. If the dissipation strength is in- 
creased beyond the "equilibrium value" , the filaments be- 
come denser and denser, and clumps in filaments may be 
formed. If finally the dissipation dominates completely the 
heating process, hot clumps, collecting almost all the mat- 
ter of the simulation zone are formed out of the filaments. 
Figs, show the change of the structure morphology for 
three different radial friction coefficients, i.e., dissipation 
strength. The friction coefficients are, Cx = 70x lO^^^ r^^l, 

= 140 X 10-3 r-i and = 210 x 10-^ r^-i. To express 
the relative strength we call the corresponding dissipations 
"weak" , "middle" and "strong" . All three simulations re- 
veal a fast fragmentation and structure formation. After 
one rotation around the galaxy center the characteristic 
striations appear already. The "weak" dissipation leads 
to a statistical equilibrium of the structure. This statisti- 
cal equilibrium establishes after about 5 rotations and is 
maintained for the rest of the simulation. It has a persis- 
tent pattern, formed by transient structures. Contrary to 
the "strong" dissipation, where the structure evolution is 
dominated by dissipation, i.e., the dissipated energy can 
not be compensated by the heating mechanism, where no 
statistical equilibrium arises 

In order to characterize more precisely the resulting 
structures we compute with Fourier transforms the 2D au- 
tocorrelation function of the matter distributions shown 
in Figs. |]^. The result is shown in Figs. |]|. The fig- 
ures reveal clearly the different morphology resulting from 
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Fig. 11. The velocity dispersions ax (a) and ay (o) 
resulting from the 2D simulations (model 1) with the 
"weak", the "middle" and the "strong" dissipation, re- 
spectively. 

the simulation with the "weak" and the "middle" dissipa- 
tion. Whereas the autocorrelation function reveals stria- 
tions with a characteristic inclination for the simulation 
with a "weak" dissipation, these striations disappear for 
the "middle" and the "strong" dissipation. 

Because we deal with self-gravitating systems which 
have negative specific heat for a certain energy range, en- 
ergy dissipation does not mean necessarily a system cool- 
ing. Indeed, the "weak" , "middle" and the "strong" dissi- 
pation cool the corresponding system only during the first 
rotation. Then the systems are heated up. After some ro- 
tations heating and energy dissipation are balanced out 
and the velocity dispersions ax and ay reach a more or 
less stable level (see Fig. [ll] ). 

It is interesting to note that ax > ay always holds for 
the "weak" and the "middle" dissipation strength. The 
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Fig. 5. The evolution of the particle positions, resulting from model 1 with a "weak" dissipation. The number of 
rotations of the shearing box around the galaxy center is indicated at the top of each panel. Shown is each second 



particle. Full resolution figures available at tittp://obswww.unige.ch/Preprints/dyn_art.htnil 
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Fig. 6. The evolution of the particle positions, resulting from model 1 with a "middle" dissipation. 



same holds also for the "strong" cooling during the first 
six rotations, then this ordering is obviously destroyed by 
the formation of hot clumps. 

In order to characterize the structure of the simulation 
terminal phase we determine the structure dimension D 
and the index S. The longer term evolution of the struc- 



tures may be superimposed by fluctuations^] on time-scales 
of the order of ~ 1/2 Trot, where Trot is the time for a ro- 
tation around the galactic center. In order to eliminate 
these fluctuations we indicate in this paper mean values 



^ To obtain an idea about these fluctuations see the evolution 
of the Schwarzschild velocity ellipsoid in Fig. |l^ and Fig. |2^. 
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Fig. 7. The evolution of the particle positions, resulting from model 1 with a "strong" dissipation. 




Fig. 8. The autocorrelation function of the structures shown in Fig. ^, resulting from the simulation with the "weak" 
dissipation. 



of the structure dimension D and the index S, determined 
during the last two rotations. In figures showing D ot 6 
we indicate in addition la error bars. 

The dimension D and the index 6 resulting form the 
structures shown in Fig. |[]^ are plotted in Fig. [ij and 
Fig. |l^, respectively, as a function of the scale R. The 
vertical lines at the left of the curves are the la error bars. 
Contrary, to the structure dimension D the error bars of 



the index 5 can vary up to a factor 7. Thus we indicate in 
Fig. |l^ in addition the position and the size of the largest 
and the smallest error bars. 

The stronger the dissipation is, the more filamentary 
resp. clumpy are the resulting structures and consequently 
the lower is the structure dimension. A comparison of the 
structure dimensions D with the indices of the velocity 
dispersion-size relation S shows that {6{R)) increases with 
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Fig. 9. The autocorrelation function of the structures shown in Fig. ^, resulting from the simulation with the "middle" 
dissipation. 




Fig. 10. The autocorrelation function of the structures shown in Fig. resulting from the simulation with the 
"strong" dissipation. 



decreasing {D{R)), where I < R < L (resolution: I — Ix — 
ly, box size in the plane: L = Lx = Ly). 

For the strong dissipation the mass-size relation can be 
approximated by a power-law for a scale range of roughly 
1 dex, but the error bars are relatively large and the scale 
range is too small to call the corresponding structure scale- 
free or fractal. 



As long as the structures are not completely dom- 
inated by hot clumps (strong dissipation) the velocity 
dispersion-size relation may be approximated by a power- 
law. However, also here the mean error bars are quite large 
with respect to the value of 6, especially for the weak dissi- 
pation, where the resulting 5-value would also be compat- 
ible with an uncorrelated velocity dispersion-size relation. 
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Fig. 14. The evolution of the particle positions, resulting from a simulation performed with model 2. The softening 
length is, e = 0.02 Acrit- Shown is each second particle. 




log(R) IX^J 

Fig. 12. The structure dimension D as a function of the 
scale R. The corresponding structures result from the 2D 
simulations (model 1) with the "weak", the "middle" and 
the "strong" dissipation, respectively. The error bars at 
the left of each curve indicate the mean Itr errors. 



For the middle dissipation there is a little more ev- 
idence for a power-law relation over roughly 1 dex. 
However, also if there is a power-law relation the value 
of 6 would be far away from the index of Larson's law 
measured in molecular clouds (0.3 < < 0.5). 

The softening length in model 1 is quite large and cor- 
responds to those of the TK model. In model 2 (Table 3 
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Fig. 13. The index 5 as a function of the scale R. The 
corresponding structures result from the 2D simulations 
(model 1) with the "weak", the "middle" and the "strong" 
dissipation, respectively. At the left of each curve the la 
mean error bars are indicated. Furthermore, the location 
and the size of the maximal and the minimal error bars 
are indicated. 



and 4) we reduce the softening length s, but we pay atten- 
tion that e > Ix = ly is always valid, i.e., that the softening 
length is always larger than the cell size of the simulation 
mesh. As expected the general tendency is that a smaller 
softening yields a stronger clustering and thus a smaller 
structure dimension. The structures resulting from sim- 
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Fig. 15. The structure dimension D as a function of the 
scale R. The corresponding structures result from a 2D 
simulation performed with model 2. The softening length 
is, e = 0.02 Acrit- The dynamical range of the simulation 
is 2.5 dex. 



ulations with a small softening length (e < 0.1) become 
relatively fast very clumpy (after 2-3 rotations). In con- 
trast to simulations with a strong dissipation, where also 
clumps are formed, the number of clumps remains from 
a certain moment on nearly constant, i.e., it does not de- 
crease due to mergers. This is shown in Fig. |lj. The struc- 
ture dimension for this simulation is shown in Fig. |l5|. 

5.2. 3D Simulations 
5.2.1. Isotropic Kernel 

We extend the models to 3 dimensions and carry on using 
an isotropic particle potential (model 3 and 4). 

In models using a particle-mesh method the number 
of particles is determined by the number of mesh-cells 
and vice versa. Due to computational limits and in order 
to do a reasonably sized parameter study on the avail- 
able machines we have to limit the number of particles to 
N « 130000. Thus it is not possible to resolve the sys- 
tem vertically with a softening length e k, — ly. The 
models using an isotropic potential reproduce thus only 
2D dynamics in a 3D space. As expected we find there- 
fore the same parameter dependence as in model 1 and 
2, respectively. Compared with previous models we carry 
out here also simulations for a slightly larger softening 
length. These simulations reveal for a limited scale range 
approximately a velocity dispersion-size power-law rela- 
tion. Fig. |l^ shows the structures resulting from 2 simula- 
tions with e = 0.3 Acrit- The indices 5 resulting from these 
2 simulations are shown in Fig. |l^. 




Fig. 16. The particle positions after 8 and 10 rotations, 
resulting from model 3 and model 4. For model 4 we only 
show each forth particle. Upper panels (model 3): Cx — 
140 X 10-3 T-i, e = 0.3 Acrit, N = 32760. Lower panels 
(model 4): = 110 X 10-^ r-,i, e = 0.3 Acrit, N = 
131040. 
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Fig. 17. The index 5 as a, function of the scale R. The cor- 
responding structures result from simulations performed 
with model 3 and 4. The index resulting from model 3 
shows a power-law relation with 5^ ~ 0.1. 



5.2.2. Anisotropic Kernel 



We use now the anisotropic kernel described in Sect. 3.3.2 



Thus the softening length e is here no longer a free pa- 
rameter. The softening of the gravitational potential is 
given by the resolution of the simulation mesh. With the 
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Fig. 18. The velocity dispersions in the plane a^y and 
vertically to it, cr^, as a function of the vertical friction 
coefficient Cz for model 7, 14 and 18 (upper panel, middle 
panel). The lower panel shows the disk scale hight zq. The 
parameter values are mean values ascertained during the 
last two rotations. As soon as the dissipation leads to the 
formation of clumps the disk is heated up. The structure 
of model 14, eg., becomes clumpy for = 50 x 10^'^ Tosc- 



anisotropic kernel we can resolve the system vertically, so 
that there is a vertical dynamical-range of 1.8 dex for all 
models with anisotropic kernel. Since the third dimension 
is now resolved we explore the structure in dependence 
on the friction coefficient Cz (model 7) and the vertical 
frequency v (model 8). 

Vertical Friction Coefficient Cz- We start with simulations, 
where Cx — Cz- However, these simulations show that 
with such a dissipation it is not possible to maintain simul- 
taneously strong density fluctuations and the disk thick- 
ness. That is, either the density fluctuations are main- 
tained and the disk scale height tends towards zero, or 
the disk scale height is maintained nearly constant and 
the density fluctuations smear out. Therefore we choose 
for all further simulations Cx > Cz, i.e., tcooi.x < ^cooi.y 
These results justify a posteriori the choice of two fric- 
tion coefficients. The velocity dispersions in the plane 
Cxy = (c^ + ^yY^'^ are mainly controlled via Cx, whereas 
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Fig. 19. The disk scale height zq and the vertical velocity 
dispersion CTz as a function of the vertical frequency ly 
deduced from model 8 and model 15, respectively, zq and 
az are mean values calculated during the last two galaxy 
rotations. 



(Jz and with it the disc scale height zq is principally driven 
by Cz- However, as soon as structures are formed with 
sizes comparable or smaller than the disk thickness, the 
dynamics in the plane and vertically to it are no longer 
independent. Thus the effect of Cz on the structure de- 
pends also on Cx- The general effect of Cz on the self- 
gravitating disk is the following: As long as the structure 
remains filamentary an increase of Cz, diminishes zq and 
(Tz- The solid curve in Fig. resulting from a simula- 
tion performed with model 7, represents such a behavior. 
The effect of Cz is also studied in model 14 and model 18. 
There the structure changes from filamentary to clumpy 
due to an increase of Cz- As a consequence these systems 
may be heated up by further dissipation (see Fig. [l8| ). 

Vertical Frequency v. In model 8 and 15 we study the 
effects of the vertical frequency v on structure and dy- 
namics. The vertical frequency determines the strength 
of the backward force due to a displacement from the 
galaxy plane. The backward force stems from the external 
galactic potential. Thus an increase of v binds the particle 
stronger to the disk and diminishes zq (see Fig. p^ ). 

In all 3D models the mean particle- velocity vertical to 
the plane (vz) is not exactly zero, but oscillates with an 
amplitude of « 0.1 kms^^ and a frequency equal to the 
vertical frequency i^. 

Concerning the effect of v on the structure, there are 
frequencies producing a clumpy structure, whereas some 
higher or lower frequencies produce a more filamentary 
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Fig. 20. The minimal structure dimension -Dmin as a 
function of u. The structure dimension is a mean value 
determined during the last two galaxy rotations. The cor- 
responding structures result from simulations performed 
with model 15. 



structure. To show this we determined the minimal struc- 
ture dimension, 



{min[D(i?,)] -.IkRkL} 



(33) 



where L = = Ly (box size in the plane) and I ~ Ix = ly 
(resolution) . The minimal structure dimension determines 
how strong the structures differ from a homogeneous mat- 
ter distribution, i.e., the lower Z^min the more filamen- 
tary resp. clumpy the structure. We calculate i^min for 
the structures resulting from simulations with different i^. 
The result is shown in Fig. The two simulations with 
^min ~ 1-7 have a more clumpy structure than the ones 
resulting from the other simulations, which arc more fila- 
mentary. 



Particle Number A^. With model 11 we determine the 
structure dimension and the disk scale height as a func- 
tion of the particle number N. We only alter the particle 
number. The mesh resolution and the friction coefficients 
Cx and Cy are kept constant. The 2D simulations of TK 
shown that a higher particle number leads to a stronger 
dissipation. Thus we expect a decrease in the disk scale 
height and a smaller structure dimension for higher N. 
We find a clear decrease of the disk scale height for an 
increasing particle number, but there is no clear trend of 
the structure dimension. This is because the disk is heated 
up in the plane due to the decreasing disk thickness. 

We find that the effect on the structure due to a change 
of the particle number can be compensated with an ap- 
propriate choice of the dissipation strength, i.e., if we use 
a weaker dissipation for an increased particle number, the 
statistical properties of the resulting structures remain un- 
changed. 

Radia l Fr iction Coefficient Cx / Large Simulation Zone. In 
Sect. 5J we explored for the 2D models how the structure 
formation depends on the radial friction coefficient Cx- 
Here we study the connection between structure and radial 
friction in 3D. Considered are structures resulting from 
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Fig. 24. The velocity dispersion components ax, Cy and 
(Tz as a function of time t (indicated in galaxy rotation 
units) . The dispersions result from a simulation performed 
with model 10 and "weak" dissipation. 
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Fig. 25. The vertical dispersion cr^ as a function of the 
dispersion in the plane axy (model 10, "weak" dissipation). 
This plot shows clearly how the system attains a stable 
state. The star and the circle denote the starting and the 
ending point respectively. 



models with large simulation zone and anisotropic kernel 
(model 5, 6, 9 and 10). 

Concerning the structure in the plane we find qualita- 
tively the same dependence as in the 2D simulations. That 
is, the initially formed structures smear out if we don't 
dissipate energy. However we still find correlations in the 
matter distribution after ten galaxy rotations. For an in- 
creasing radial dissipation the structures become denser 
and denser until finally hot clumps are formed out of fil- 
aments. Fig. ^ and |2^ show, as an example, the struc- 
tures resulting from a simulation performed with model 
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Fig. 21. The evolution of the particle positions seen from above the galaxy plane. The structures result from a 
simulation of model 10. The friction coefficient is = 64 x 10~^ ("weak" dissipation). The number of rotations 
of the shearing box around the galaxy center is indicated at the top of each panel. Shown is each fourth particle. 




Fig. 22. The particle positions inside the slice, — O.OSAcrit < U < O.OSAcrit, seen along the direction of orbital motion 
(model 10, "weak" dissipation). 
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Fig. 23. The autocorrelation function of the particle distribution, presented in Fig. [2^ and ^ (model 10, "weak" 
dissipation). Upper panels: Autocorrelation function in the x-y-plane. Lower panels: Autocorrelation function in the 
z-x-plane. 



10. The structures reach in this simulation very fast (after 
2 rotations) a statistical equilibrium. The autocorrelation 
functions revealing the underlying characteristics of these 
structures are shown in Fig. Fig. ^ and show the 
evolution of the velocity dispersions. These Figures con- 
firm that a statistical equilibrium is attained after w 2 
rotations around the galaxy center. 

The filaments resulting from these simulations {Cx — 
64 X 10~^ Tosc) are very dense and they are first signs 
of clump formation. A slight increase of the dissipation 
would thus turn the structure from filamentary to clumpy. 
Indeed, the structures in Fig. ^resulting from simulations 
performed with the same model but with slightly higher 
radial dissipations Cx = 66, 68, 70 x 10~^ Tqsc, are already 
clumpy. For convenience we call the relative strength of 
the radial dissipations used in model 10 "weak" , "mid- 
dle" , "strong" , and "very strong" . In Fig. ^ the structure 
dimension D is shown for the four different dissipation 
strengths. The structure dimensions are not constant over 
the corresponding dynamical range and are thus not frac- 
tal. However, the structure dimension resulting from the 
simulation with the "strong" dissipation has a dimension 

1.5 < D < 1.8 over the whole dynamical range in the 
plane and remains smaller than 2 also on scales where the 
disk thickness becomes important. Furthermore, the struc- 
ture dimension has at 0.06 Acrit the same dimension as at 

1.6 Acrit with a minimum in between. This is qualitatively 
a different behavior than those of the initial state. The 
slight increase of the structure dimension at the left and 
at the right may be induced by the small scale (resolution) 
and large scale cutoff (periodicity), respectively. This be- 



havior is quite general for our simulations and it may be 
that a larger dynamical range produce a more constant 
structure dimensions, D{R) « const.. 

We do not find a velocity dispersion-size relation sim- 
ilar to a power-law for the simulations performed with 
model 10. However we do find some hints for such a re- 
lation in the structures resulting from model 6. Fig. [2^ 
shows the velocity dispersion-size relation for two differ- 
ent simulations of model 6. The simulation with the weak 
dissipation produces filamentary structures whereas the 
simulation with the strong dissipation produces filaments 
and clumps, thus S varies stronger and the error bars are 
larger for these structures. However, both simulations pro- 
duce a velocity dispersion-size relation which may be ap- 
proximated to first order by a power-law. Of course, the 
error bars are too large and the scale range is too small 
to call the corresponding velocity dispersion-size relations 
scale free. The resolution in model 6 is larger than those in 
model 10. Why this softer gravitation seems to reproduce 
better a power-law velocity dispersion-size relation is at 
the moment unclear. 

Analog to the 2D case we find also in 3D simulations 
a systematic ordering of the velocity dispersion compo- 
nents. The ordering, ax > cFy > az, holds for quite a 
large range of dissipation strength. This is shown by sim- 
ulations performed with model 10. The velocity disper- 
sion components resulting from the simulation with the 
"weak" dissipation strength is shown in Fig. |4| and those 
with the "middle" , "strong" and "very strong" dissipation 
strengths are shown in Fig. |2^. Even if we don't dissipate 
energy this ordering is still observed after 10 galaxy ro- 
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Fig. 26. The particle positions after 8 and 10 galaxy 
rotations. The structures result from simulations per- 
formed with model 10. Upper panels: "middle" dissipa- 
tion. Middle panels: "strong" dissipation. Lower panels: 
"very strong" dissipation. Shown is each forth particle. 



tations. Only strongest clump formation can destroy the 
sys tematic ordering. From the definition of a; and y in Sect. 
3.1 it follows ax = CTfi and <7y = a^. To sum up, one can 
therefore say, that as long as our model produces struc- 
tures resembling the lumpy matter distribution in spirals, 
it also produces an ordering of the velocity dispersion com- 
ponents (aji > (Jif, > (Jz) corresponding to those of the so- 
lar neighborhood, or of TV-body simulations of spiral disks 
(Pfenniger & Friedh 1991). 



Radial Friction Coefficient Cx / Small Simulation Zone. In 

order to increase the resolution in the plane and to ap- 
proach the vertical resolution we decrease the size of the 
simulation box. The structures resulting from these sim- 
ulations differ from those in the large simulation box. We 
still find filaments for a "weak" dissipation and clumps 
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Fig. 27. The structure dimension D as a function of the 
scale R. The corresponding structures result from 3D sim- 
ulations of model 10. The structure dimension is indi- 
cated for the simulation with the "weak" (see Fig. ^l]) , 
the "middle" , "strong" , and the "very strong" dissipation 
(see Fig. H). 




Fig. 28. The index (5 as a function of the scale R. The 
corresponding structures result from 2 simulations per- 
formed with model 6. One simulation is performed with 
a "weak" dissipation and produce a filamentary structure 
whereas the other simulation is performed with a "strong" 
dissipation, producing a clumpy structure. 



for a "strong" dissipation (see Fig. |3^). However, filaments 
and clumps, respectively, appear less numerous. From this 
point of view the structures in the small simulation box 
are not a fractal continuation of those in the large simu- 
lation box. However, this does not exclude the possibility 
that a simulation with a dynamical range incorporating 
the scale of the large and the small simulation box would 
produce scaling laws over the whole dynamical range. 
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Fig. 30. The evolution of the particle positions seen from above the galaxy plane. Shown is each forth particle. The 
structures result from simulations performed with model 16. We show 3 simulations with different dissipation strength. 
Upper panels: "weak" dissipation. Middle panels: "middle" dissipation. Lower panels: "weak" dissipation. The number 
of rotations of the shearing box around the galaxy center is indicated at the top of each panel. 



Shear. The inhomogeneous structures appearing in the 
shearing-box simulations can only be maintained when the 
dissipated energy is compensated by an energy injection. 
The source of this energy is the shear motion. If the shear 
is reduced, the dissipated energy can no longer be bal- 
anced by the energy injection and the system collapses. 
A rotation curve, increasing with the square root of the 
radius {v (x -y/r) reduces the shear-flow, y — —2Aox, by a 
factor two. With model 19 we perform simulations for dif- 
ferent rotation curves. That is, for the usual flat curve with 
Aq = 0.5 JIq, k = \/2 ^^0 SiJid for a curve increasing with 
the square root of the radius, Aq — 0.25 fioi k = V^flo- 
Such a choice of the rotation curve doesn't reflect a realis- 
tic case but serves to explore the influence of the shear on 
the structure formation. However, a change of the rotation 
curve alters also the epicyclic frequency and consequently 
the Schwarzschild velocity ellipsoid as well as the critical 



wavelength Adit • Because we only want to examine the ef- 
fect of the shear, simulations with different rotation curves 
are carried out with the same initial velocity dispersion. 
Moreover, the distances are indicated in units of kpc. The 
effect of a decreased shear is revealed in Fig. |3^. The same 
friction coefficients as in the simulations with the flat rota- 
tion curve lead for the slightly increasing [v oc ^/r) curve 
to a collapsed system. 

5.2.3. Minimal Structure Dimension and Maximal 
Index 5 

In this paper we study mainly the structure dimension in 
dependence of the radial friction coefficient Cx, because 
this parameter determines principally the degree of struc- 
ture inhomogeneity. Indeed, the parameters and serve 
mainly to prevent the disk from vertical dissolution and a 
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Fig. 31. Evolution of the matter distribution, seen from above the galaxy plane. The structures result from simulations 
with different rotation curves performed with model 19. Upper panels: flat rotation curve {v = const.). Lower panels: 
increasing rotation curve (v ex ^/r). The number of rotations around the galaxy center is indicated at the top of each 
panel. Shown is each forth particle. 



change of the differential rotation or the particle number 
can be balanced by an appropriate choice of the dissipa- 
tion strength. 

The structure characteristics in dependence of Cx are 
studied for different simulation boxes and resolutions. 
In order to compare the structures and the dynamics 
as a function of Cx resulting from the different models 
we calculate the minimal structure dimension and the 
maximal index S of the velocity dispersion-size relation. 
The minimal structure dimension is defined in Eq. (p3|). 
Correspondingly, the maximal index is, 



{max[(5(i?)] ■.1<R<L} 



(34) 



The results are summarized in Fig. The general trend 
is that the structures become denser and more clumpy for 
higher radial dissipation, leading to a smaller structure di- 
mension and to higher velocity differences on the different 
scales of the system, i.e., to higher Jmax- 

6. Discussion 

With respect to the models with the large simulation zone 
(model 1-11), the models with the small simulation zone 
(model 12-19) produce less fragmentation. That is, only a 
few weak filaments appear in simulations of these mod- 
els as long as the dissipation is weak and an increase 
of the dissipation strength can not produce a more frag- 
mented or filamentary structure, because clumps appear 
very quickly. These clumps become rapidly unphysically 



hot and collect almost all the matter of the simulation 
zone. The clumps may thus hinder the formation of a more 
fragmented structure in these models. 

The formation of non-transient collapsed clumps, 
which have the tendency to attract more and more matter 
were also found in other numerical studies using dissipa- 
tive or sticky particles in order to model the ISM (Huber 
^00 1| , Semehn & Combes [2000| ). The appearance of these 
clumps may thus be a generic problem of gravitational 
clustering simulations with dissipative particles. 

In the following we discuss some aspects related with 
the formation of these non-transient and in our opinion 
unphysical hot clumps, which may hinder the evolution 
towards a structure, being hierarchical over a more ex- 
tended scale-range. We discuss numerical, dynamical, and 
physical aspects of the problem. 

6.1. Numerical Problems 

The size of dissipative particles, given by the softening 
length, is much larger than the smallest clumps in the 
interstellar medium. Thus the nearly homogeneous mass 
distribution (smear out) over the softening length is not 
justified and it is not a priori clear that the inhomoge- 
neous structures below the resolution scale do not affect 
the structures on larger scales. 

Observations of the ISM reveal filamentary structures 
down to the smallest resolvable scales, thus it is unclear 
down to what scales the highly inhomogeneous structures 
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Fig. 29. The velocity dispersion components Ox, <Jy and 
(Tz as a function of time. The dispersions result from the 
simulation with the "middle" , "strong" and "very strong" 
dissipation strength, respectively (model 10). The matter 
distribution of the corresponding simulations are shown in 
Fig. |§. 



continues. Moreover, a large scale range has probably to 
be taken into account in order to reproduce the observed 
structures. Thus it is an open question whether in the fol- 
lowing years a resolution will be reached such that basic 
clumps can indeed be represented by dissipative particles. 
Taking into account the ubiquitous trend of gravitation- 
ally unstable media to produce sheets and filaments not 
only in the ISM but also in cosmological simulations, it 
might be that the particle model of the basic mass unit, 
as a rigid body conserving mass, is not adequate. 

6.2. Dynamical aspects 

One could argue that when the "mean free path" of cloud 
clumps is larger than their size and a particle description 
seems to be justified, physical cloud clumps collide and dis- 
solve because they contain internal degrees of freedom due 
to the smaller subclumps moving inside them. Thus dis- 
sipative particle simulations should incorporate collisions 
with mass exchange (Pfenniger 1994 ). This is particularly 
relevant if the clouds collide with supersonic speed, that 
is, with velocities larger than their internal velocities. 
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Upper panel: The minimal structure dimen- 
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Fig. 32 

sion D 

Cx ■ Lower panel: The maximal index Jmax of the velocity 
dispersion-size relation as a function of the radial dissipa- 
tion strength. Dmin and ^max are mean values calculated 
during the last two galaxy rotations. 



Inherent dynamical properties of gravitational unsta- 
ble media may cause further problems. Let us discuss two 
of them: 

1.) Typically clumps are subject to an anisotropic grav- 
itational contraction that alter its morphological charac- 
teristics, i.e., clumps may become pancakes which become 
filaments. This was shown by Zeldovich (1970) and nu- 
merically confirmed by Kuhlman et al. (1996, hereafter 
KU). In their numerical experiments KU replaced parti- 
cles in dense regions by clouds, made up of 2^ resp. 2^ 



particles. They found that only a small fraction collapses 
along all three directions and forms dense clumps. A fur- 
ther subdivision of the particles would probably lead to 
the same result. If the transformation from clumps to pan- 
cakes and to filaments on small scales is important for the 
appearance of the global structure, it may be problem- 
atic to model gravitational unstable media with particles, 
because a higher particle number would not solve the in- 
herent problem of anisotropic clustering. 
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2. ) A further problem when simulating gravitational clus- 
tering is related to the exponential propagation of two- 
body relaxation in hierarchical A^-body systems. If the 
mass distribution of the considered system is self-similar 
over a range of scales, at each scale the effective bodies 
can be viewed as the corresponding clumps. The hierar- 
chical structure acts as a strong two-body relaxation am- 
plifier since at each scale the effective number of bodies 
is strongly reduced. If this effective number is 0(10) the 
relaxation time at each level is of the order of the crossing 
time, so two-body relaxation is a major driver of evolution 
throughout the scales. 

3. ) Let us describe in more detail the related problem of 
error propagation in a hierarchical system. First we show 
that in a gravitational unstable medium developing a hier- 
archical structure, matter on smallest scales evolves faster. 

A hierarchical mass distribution satisfies the scaling 
relation between two levels L and I: 

D 

(35) 



Ml 
Ml 



Rl 
Ri 



where D is the mass scaling exponent, which would be the 
fractal dimension in a self-similar hierarchical system. D 
is restricted to stay in the interval [0 — 3] since mass is 
positive and space filling can not exceed the third power 
of scale. Then the density scales as 



PL_ 

Pi 



Rl 
Rl 



Z?-3 



and consequently the crossing time Tdyn oc [Gp) scales 
as 

^ = (t) . (37) 

Tdyn,; \^l J 

So in a hierarchical model the dynamical time (or crossing, 
or free-fall time) always decreases a t sma ller scales since 
< £» < 3 (Pfenniger & Combes |l994[) . Thus the low 
scale structures evolve faster. 

Self-gravitating A^-body systems are chaotic and 
neighboring trajectories diverge exponentially (Miller 
1964). This means also that any error propagates expo- 
nentially: 



Ax oc Aa;oe 



At 



(38) 



where Aa;o is a small initial error at time i = 0, and Ax 
is the error at time t. Now, the degree of chaos and thus 
the error evolution is determined through the maximum 
Liapunov exponent A, which for small A^-body systems 
is approximately inversely proportional to the dynamical 
time, A oc 1/Tdy„ (Miller |1994| ). 

With this estimate, let us see how errors are amplified 
through the two adjacent levels. If Ax; is in fact the initial 
error at the lower level I, we get after one crossing time at 
level L, Tdynx: 

(3-D)/2\ 



AXL 
Ax; 



exp 



( "^dynj 
V Tdyn,; 



exp 




The error amplification becomes rapidly huge as soon as 
D <i. Say, if = 2 and D = 1.5, then ^ = 5.38. 
Across n levels the error ratio at the highest level goes as 
5.38" and becomes much larger than the scale ratio, 2". (It 
is easy to show that for D < 2.264 the error amplification 
is always larger that the scale ratio, while for > 2.72 
and D > 2.265 one can find hierarchical systems for which 
this problem does not occur). 

Thus we find an inherent dynamical problem, which 
can not be solved by using a higher resolution, because 
the increase of resolution exponentially increases the er- 
rors through the scales. Possibly such hierarchical systems 
might be dominated by numerical errors in simulations, 
and by small scale physics in real systems. 



6.3. Additional physics 

Small scale physics, supporting the dissolving process of 
dense clumps, is not taken into account in our simple 
model. However, stellar winds, supernovae, jets and out- 
fiows may be important in the overall mass transport 
across the scales. These processes may ensure a cyclic 
matter flow by giving matter back to larger scales, which 
would condense back via gas cooling. Such a matter-flow 
may be crucial for sustaining the transient nature of hier- 
archical clumps for extended time. 



(36) 

7. Conclusions 



(39) 



The structure resulting from the local simulations of self- 
gravitating disks can be homogeneous, filamentary or 
clumpy depending on the relative strengths of the compet- 
ing gravitational and dissipation processes. As long as the 
structure is mainly filamentary self-gravitation and dis- 
sipation ensure a statistical equilibrium, where repeated 
transient patterns are formed. If the dissipative processes 
begin to dominate the evolution, the structures turn from 
filamentary to clumpy. During the subsequent evolution 
the clumps become hotter and more massive. 

In general, clumpy structures do not evolve towards a 
statistical equilibrium. However model 2 does show that it 
is also possible to establish a persistent pattern of clumps. 
These clumpy structures may show signs of a power-low 
mass-size relation. Some of our 2D as well as 3D simu- 
lations suggest also a power-law velocity dispersion-size 
relation. 

However the scale range of the simulations appears too 
small to draw final conclusions about a fractal structure 
and an extension of Larson's law beyond the size of molec- 
ular clouds. We can suggest a few reasons causing the dis- 
crepancy: 

1. ) The numerical resolution should extend over several 
more decades of scale before the fragmentation stops to 
be dominated by finite scale range boundary effects. 

2. ) A fundamental law is associated with the rigid point 
particle representation of mass which, by forcing a partic- 
ular lowest scale boundary condition, would prevent the 
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bottom-up building of scaling relations to match the ob- 
served ones. Also the propagation of errors may be super- 
exponential in a hierarchical organized medium because 
the error evolution at largest scales aren't determined 
by the dynamical time of these scales, but by the much 
smaller dynamical time of the smallest scales. Systems 
with dimension lower than about 2.2 are particularly con- 
cerned. This point is relevant for particle simulations of 
gravitationally unstable systems, such as cosmological and 
thin disk {D < 2) simulations. 

3.) A key physical ingredient, such as mass cycling through 
the scales due to star formation, stellar activity, and gas 
cooling could be essential to sustain a fractal state of the 
ISM. 

Finally, it is interesting to note that the anisotropy of 
the velocity-dispersion ellipsoid, resulting from our simu- 
lations, has systematically the same ordering ((T/? > > 
Uz) and relative amplitudes as observed in the Galaxy 
and in A'^-body simulations of spirals. Since the models 
are deliberately a simplified representation of reality, we 
learn from this that this ordering may be due to a very 
general property of galactic disks, to be substantially self- 
gravitating in 2, and to rotate differentially with a similar 
shear rate set by a constant rotation curve. 



Appendix A: Pseudo code of the shearing box 
program 

• Initialization of the particle positions and velocities. 

• Calculation of the canonical momenta. 

• The inclination of the meshes is 
flb = [(— rioO^od (Ly/Lx)] and 
a/ = Of, -I- {Ly/Lx), respectively. 

Thus there is a finite number of possible inclinations, 
n = T/At, where T = Ly/{Lx^o) and At is the time- 
step. 



The particle positions of the Cartesian coordi- 
nates are saved {yc = y). 

Transformation of the particle position: Cartesian 
coordinates Forward mesh {y ~ y — ax). 



else 



a = a,. = af — Ly/Lx- The inclination a corre- 
sponds to the backward mesh. 

y = yc- The particle positions are again those in 
Cartesian coordinates. 

y — y— ax. Transformation of the particle position: 
Cartesian coordinates — s- Backward mesh. 
The missing points in the V-fC(af,)-matrix are de- 
duced from VAr(a/), making use of their symme- 
try. 



endif 

• Transformation of VA'(a) 
{VK{a) VK{a)) 
mod L„ 



into Fourier space. 



y = y mod-Lj,. All the particles should lie inside the 
local box of the forward and the backward mesh, 
respectively. 

Mass to mesh assignment and transformation of the 
density p into Fourier space (p ^ p). 

of the potential differential (V$ = 



), where ^ denotes the inverse Fourier 



Calculation 
transform. 

Calculation of the forces acting on the nodes of the 
affined meshes. 



F = a— - — 

fill riT 



9y 
-(1 

a<s> 



dx 



dx 



• Force interpolation in accordance with the CIC- 
method. 

• Transformation of the forces. Affine mesh 
Cartesian coordinate system. 

enddo (5 = 1,2) 

• The particle positions are again those in Cartesian 
coordinates (y = yc) 

• Force weighting 

• Calculation of the new canonical momenta and parti- 
cle positions by means of the implicit canonic al fini te 
difference approximation (Pfenniger & Friedli 1993| ). 

if (z-coordinate of a particle lies outside the local box.) 

• The particle is considered as escaped, 
do g = 1, 2 (Do for forward and backward mesh, respectively) 

elseif {x- and/or y-coordinate lies outside the local 
box.) 

a — af ^ [(—Q,ot)'mod{Ly/Lx)] + Ly/L^- The 
inclination a corresponds to the forward mesh. 



do i=l,n 

• The derivations of the Kernel \^K{af) for the differ- 
ent possible inclinations a/ of the forward mesh are 
calculated. K in an xuyX n^-matrix, where UxnyUz 
is the number of cells in the local simulation box. 

• The points: {nx/2, ny/2, 1 : Hz) of VK{ab) are calcu- 
lated for the backward mesh. The other Kernel points 
can be deduced from 'VK{af) of the forward mesh, 
making use of their symmetry. 

enddo 

do i— l,nt (time propagation loop) 



if (g=l) 
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• Reentrance at the opposite side with appropriate 
canononical momenta (see Eqs. (^). 

endif 

if (storage condition) 

• Calculation of the particle velocities. 

• Storing of the particle positions and velocities to the 
disk. 

endif 

enddo (time propagation loop) 
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